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SUMMARY 


This work outlines an upwind, second order accuracy, 
coupled, conservative numerical scheme for solving a two dimen- 
sional laminar or turbulent flow field. Mean vorticity, w, and 
mean stream function, <p , are used as the mean flow dependent 
variables. The turbulent kinetic energy k and the turbulen-c 
energy decay rate, e, are used to define the turbulence state. 

The rate of convergence of the coupled, conservative, mean field 
system as well as the turbulence state system k~e , is twice 
that realized when solving these equations separately. Although 
the turbulence boundary conditions have a non-regular variation 
near a solid wall, the turbulence model is reduced exactly to this 
variation, keeping the conservative features. 

The axisymraetric mixing of two confined jets with an internal 

heat source is considered with this numerical scheme. The inlet 

boundary conditions have a limited effect if they are applied 

far upstream of the end of the inner cylinder (" 5 radii) . The 
* 

fully developed flov; conditions which apply at the exit section, 
should be located far downstream to prevent an effective change 
in the real flow Reynolds number. The parabolic conditions at 
the exit section are found to be very good. A laminar recirculation 
zone appears f r certain flow parameters whose length tends to 
som.e asymptotic value as the Reynolds number increases. During 
the variation of the Reynolds number, the recirculation mass 
flux approaches a maximum. Thermal radiation from the heat 
generating inner jet material is the leading mechanism of thermal 
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energy generating transport. This transport does not appear to 
depend heavily on the turbulent diffusion coefficients. The 
temperature field depends mainly on the radiative conductivity 
coefficients and on the heat generation. The latter is the only 
term that depends on the turbulence status through the species 
concentrations. The conditions for obtaining stable turbulent 
flow solutions are stated in this study, and are, almost always, 


fulfilled. 


I . INTRODUCTION 


Computational methods, as a tool to analyze ohysical phenomena 
and to solve the related engineering problems, are becoming 
even more popular today. This is mainly due to the rapid growth 
in computer capability (in terms of direct access memory and the 
improvement in the CPU time per elementary action) and improve- 
ments in the basic numerical schemes. 

This report concentrates on developing computational method 
v/ith which to obtain a physical understanding of the turbulent 
field of two coaxial jets entering an axisymmetric chamber, Evan 
the laminar field of this flov/ is quite complicated. This is 
due to the many different domains which exist in the field 
especially in the entrance region. Physically, three regions 
may be identified: the wall region, the initial region near the 

axis of symmetry and the mixing region. Advancing downstream, 
these regions change relative size v;ith the ratio of the two 
jets' mass fluxes as the main parameter. The turbulent field 
of these flows is much more complicated due to the difference 
in the effective transport coefficients and turbulence level 
from region to region. However, being aware beforehand of the 
complications and the different regions of this field, one can 
adjust the appropriate turbulence model and numerical scheme to 
treat the problem. 

The objective of the present study is to describe numerically 
the velocity, concentration and temperature fields in confined 
coaxial turbulent flows with internal energy generation. 


Physically, the flow system considered has an inner flow which is 

slower and heavier than the outer flow. The inner flow material 

undergoes nuclear fission in the chamber cavity and generates 

heat. This system is representative of a gas core nuclear reactor 

where the inner gas is fissionable and the outer flow is some 

lighter working fluid. The resultant flov/ is usually turbulent. 

In applying a turbulent model to a flow one has to compromise 

between choosing a high order of closure and paying in an 

extremely complex code for computing the phenomena, and choosing 

a very simple algebraic model and paying in a poor prediction 

of the phenomena. In this study the "two equation" model of the 
[ 1 ] 

type k-e (k is the turbulence level and e is the turbulence 

dissipation) was adopted, since, as has been stated in the 
. [ 2 ] 

literature , it is the closest-to-optimal model to predict 
two dimensional elliptic fluid flows available today. Both 
k and e are governed by convection-diffusion like equations. 

Once k and e are known at every point in the field, the effective 
transport coefficients can be calculated with models for the 
turbulent Prandtl cind Scfimidt numbers. Since the rate of energy 
generation in the core is very high, very high temperatures occur. 
Therefore thermal radiation is the dominant mechanism of energy 
transport (and the conductive transport is not negligible) . 

In such a case it is very convenient to apply the Rosseland energy 
diffusion approximation . This approximation, v/hich is valid 

only for flows with high absorptivity results in a radiative 
conductivity like term in the equations. 
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The following variables define the two dimensional axi- 
symmetric field: 

2 velocity components (u,v) 

2 rurbulent state variables (k,e) 

1 temperature variable (T) 

1 pressure variable (P) 

1 concentration variable (C) 

7 variables 

These variables are controlled by the following governing 
equations : 

1 continuity 

2 momentum 

2 turbulence model 
1 energy 
1 species 
7 equations 

In the present field, as well as in many other fields it is very 
convenient to choose the stream function (ij;) -vorticity (w) 
variables instead of the pressure (p) -velocity (u,v) variables, 
because the number of equations and unknowns is then reduced 
to six (6) . 

The nuir.erical solution is performed on a non-uniform mesh 

grid that is spread over the domain of interest. The finite 

differences for all the terms are accurate to the second order, 

[4 51 

and unlike some other schemes ' ' , no artificial viscosity 
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is inserted and the flow is so3.ved for the correct Reynolds 
number. In order to improve the stability, the equations are 
solved with as much coupling as is possible. The flow field 
equations (>|»-w) are solved in a coupled manner, as well as the 
turbulent variables (k-e) . The energy equation and the 
concentration equation are solved separately. The algebraic 
difference equations are solved iteratively by successive line 
relaxation [or generally, by the successive block line relaxation] . 
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2. FORMULATION OF THE PROBLEM. AND PHYSICAL ASSUMPTIONS 


The flow system under consideration is sketched in Fig. 1.1. 
The system consists of two concentric cylinders/ with the outer 
cylinder of diameter continuous to a distance and the 
inner cylinder of diameter < D2/ is continuous to a distance 
^1 ^ ^2' ends at the section A-A. 


j h ^ \ , 



Fig. 1.1. Schematic Fields of the Physical Problem. 

Different fluids with different mass fltixes and densities move 
through the inner cylinder (denoted with lower index 1) and 
through the annulus between the inner and the outer cylinder 
(denoted with lower index 2) . Despite the simplicity of the flow 
system/ complex physical processes take place, relating to the 
detachment of the stream from the inner cylinder trailing edge, 
contraction or expansion of the inner stream after the section A-A, 
and resulting for some cases in a recirculation region downstream 
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of the Goction A~A. Dononding on tho goometrical ratio 
the mass flux ratio 0 ,/Q ,, and the laminar Reynolds number 
ratio, a recirculation 20110 can occur either on the axis of 
synunetry or on tho outer cylindrical wall. There is also some 
weak dependence of this phenomenon on the upstream turbulent 
properties ratio. In order to describe <ho complexity of the 
physical processes occurring within the nozzle area in terms of 
the governing equations, these should be of the elliptic type. 

In the subsequent part of this report, a mathematical 
model of tlie flow in this eeometry will be presented. Lot us 
begin with the following assumptions; 

(i) The fluids '»rwi incompressible and of constant density. 

(ii) A fully developed turbulent flow exists at the 
entrance to the computational region, a distance 
before the section A-A. 

(iii) The flow properties are isothermal. 

(iv) The duct walls are impermeable. 

These assumptions will lead to a mathematical model of the flow 
field. With this model, described in the next section, it is 
possible to determine the velocity distribution in the compu- 
tational domain as well as approximate the temperature and 
concentration distributions. 
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3. MATHEmTICAL MODEL OP THE FLOW FIELD 


3 . 1 The Governlf.q Transport Equations o£ the Flow 

Steady, axially-symmetric flow of a viscous and incoippressible 

gas with constant molecular transport coefficients is described 

by the equations of continuous media motion, namely the momentixm 

equations coupled with the equation of continuity. These Navier- 

f (5 1 

stokes equations in cylindrical coordinates, r,z, are: ^ ^ 


continuity: 

3u . 1 3 
3z r 3r 


(rv) 
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(3.1) 
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(3.3) 


where u,v are the mean velocity components in the z and r 
directions, p is the density and p is the mean pressure. The 
quantity appearing in the equations (3.2) and (3.3) is the 

so-called effective viscosity, being +•’•'3 sum of the molecular 
(laminar) viscosity and the turbulent viscosity y^. It 
characterizes the turbulence status at a given point of the 
flow field: 
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(3.4) 


A turbulent flow is thus treated like a laminar flow with a 
variable viscosity coefficient, which is characterized by special 
rations to the meen flow and other turbulence quantities. This 
relation is the main result of the turbulence model. 


3 *2 The Turbulence Model 

A broad review of turbulence models is presented by Launder 
(71 

and Spalding. The various models are classified according 

to their complexity and useability, which depend on the 
particular physical hypotheses employed to describe different 
turbulent fields. On the basis of physical analysis of the 
phenomena taking place in the field of the two coaxial flows , 
a ’’two equation" turbulence model is adopted as recommended in 
.reference [2] . This model consists of tv/o dynamic equations for 
two turbulence variables: the first equation is for the turbulent 

kinetic energy (k) ; and the second equation is for the dissipation 
(e) of the turbulent energy. 

The turbulent kinetic energy is defined by the half mean 
of the sum of the square of the turbulent velocity fluctuations: 



and the turbulent energy dissipation (or decay) is defined by 


e 



9u: 
(— ) 


2 


j 


10 


OF y if 


v/here in both formulations there is a tensorial summation over 
1 <. i/ j <. 3, and is the velocity fluctuation in the i direction. 
It can be seen that k and e are positive definite quantities. 

A general form of the k-e model has been given by Launder and 
Spalding.^ The transport equations for k and e, in cylindrical 
coordinates are of the following form: 


D(u 



V 


3k ^ 

ar' 


= i_ ills + i i_ ^eff 

3s 3z' r 3r 3r' 


+ Ugff G - PC 


(3.5) 


H * n> 


/ “eff 3£l i 1 1_ ( '^eff 

3z V 3z^ r 3r a_ 3r^ 

e e 


k ^^1 ^eff ^ " ^2 


(3.6) 


where the generation term is : 


G - 2[(ai) + (-^) +(-)] + (^ + 


(3.7) 


knowing local values of k and e one may determine local values of 

f 8 1 

turbulent viscosity from the Prandtl and Kolmogorov ' ^ formula: 



The values of the constants appearing in equations (3.5), (3.6) 

and (3.8) for the turbulence model are suggested by Launder 
[2] 

and Spaldxng as follows: 
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= 1.44, C 2 = 1.92, = 0.0^ 

= 1.0, Og = 1.3 (3.9) 

These values will be reviewed later on in this report. 


3 . 3 The Energy and Species Equations 

Since the flow field in the present study is incompressible, 
the Eckert number is zero and therefore the substantial derivatives 
of pressure may be ignored in the mean energy equation. In the 
absence of body forces the following energy equation may be 
considered: 




k k' + k 


r T 3T, 


+ (}) + 


(3.10) 


where X is the total thermal diffusivity, and ({) is the dissipation 
given by 




2<k> 


2 (|) 


+ ( 


3u 

8r 


3v 

3z 


) ] 


(3.11) 


and is equal to the generation term of the turbulence energy 
given partially by eg. (3.7). Generally X is given by 

X = + X^ + X^, (3.12) 

where the first term on the right hand side of eq. (3.12) is 
the laminar contribution and the second term is the turbulent 
contribution to the thermal conductivity; X^ stands for contributions 
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from other conduction-like processes. in eq. (3.10) is the 
energy generation rate term, representing an energy source 
such as a chemical or nuclear reaction. In the present work we 
are interested especially in the effects of nuclear fission on 
the field's temperature. Regularly, the nuclear-reactive species 
enter the field through the inner pipe v/ith concentration C, and 
the fission is carried out in the core of the field. The major 
mode of transfer of the nuclear energy is by radiation processes. 
According to some well defined physical models, like the 
Rosseland diffusion approximation X^, which is a radiative 
conductivity coefficient, can have only a positive non-zero 
value, and the undimens ional source term will be of the form. 


S 


- Arc [1 - i (2 I 


4 

1 ) ] 


(3.13) 


where A is a constant and L = is the axial length of the 

chamber. The formula for X is 

r 



r 3 1 

The values of the constants were taken to be 


B 


32 

3 t3 


A :: 50.0 


(3.15) 


where t a 8 , 3 ~ 400 (®K) and the temperature T is measured in °K. 

The equation of species continuity is taken similar in form 
to that of the energy equation since both temperature and 
concentration are scalar quantities. The difference is that 
the quantity C is not generated in or lost from the system. 
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This form of the species continuity equation is an approximation 
which is only good for small variations in density. Hov/ever 
the simplification resulting from this assumption warrants its 
use in this illustrative development. The governing equation 
is then: 

k h <"!§>■ ^ !? 

where 5 is the total diffusion coefficient defined to be 

D = (3.17) 

A 

where and are the laminar and the turbulent diffusion 
coefficients . 

3 . 4 The Dimensionless Form of the Equations 

The parameters that control the flow field behavior can 
be found by casting the equations in dimensionless foirm. The 
velocities will be normalized with respect to the axial velocity 
comoonent on the center line at the entrance section, U (see 

o 

Fig .3.1). 
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The lengths are normalized with respect to the outer radius R 2 . 

Tho temperature is normalized with respect to AT = TQ-T 2 , 

where Tq is the temperature at the centerline and T 2 is the 

outer cylinder wall temperature, both at the entrance section. 

The turbulence kinetic energy and the turbulent dissipation 

2 3 

rate are normalized with respect to Uq and Uq/D 2 respectively. 

2 

(The pressure is normalized by pUq.) Therefore the flow field 
behavior is controlled by the following non-dimensional 
numbers : 

u D 

Reynolds numbers Re = -- — 


Laminar Prandtl number: 


Pr 



D 


Turbulent Prandtl number: 
Laminar Schmidt number; 


Pr = .-..li-P, 



Turbulent Schmidt number: 


Eckert number; E 


Radiation number; 



C AT 
P 







The following values for some of the above numbers will 
be adopted; 

Pr = 0.72, Sc = 0.8 

Pr^ = 1.0, Sc^ = 1.0 

E << 1 

>> 1 (20 f 100) 
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3 . 5 The Computation for an IncortipressiblG Laminar Axisynuuetric Flow 
For an incompressible, laminar flow of a Newtonian fluid with 
constant viscosity in the absence of body forces, the axisymmetric 
Navier Stokes equations and the continuity equation (3.1 - 3.3), 
retaining the time term are: 

(3.18a) 


9u 


, 9u 

1 

3P . 
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9r^ 
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3P , 

. 1 

( 2 + 
9r^ 

9^v 

. 1 

9v 

9t 

V + 

3r 

u = 

dz 

p 

9r Re 


+ — 
r 

9r 


V 


- ^) (3.18b) 

r 


(rV) + 1^ (rU) = 0 (3.18c) 

These equations are referred to as the primitive variables 
equations since they are written in terms of the basic variables 
V, u and P. If the pressure is of primary importance in a 
particular incompressible field, the equations in this form plus 
the Poisson equation for pressure instead of the continuity 
equation (3.18c) should be considered.^ ^ Experience in solving 
these equations may also be useful if one is ultimately interested 
in three dimensional flows. Solution procedures for solving the 
primitive variable equations in both two- and three-dimensional 
rectangular coordinates are discussed by many investigators . 

Since the pressure field is not of main interest in the present 
study, the following stream-function-vorticity formulation will 
be adopted. 

By eliminating the pressure term in eqs . (3.18a) and 

(3.18b) by cross differentiation, applying eq. (3.18c) to this 
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result, and then using the definition of the vorticity «: 


w 


9z 


9r 


(3.19) 


the following dynamic equation for the vorticity may be obtained: 


So) . vcj 

3 t ^Ir '^Tz"r“ 


T? 3.2 r 3r - ^2^ 


(3.20) 


This "vorticity transport" form of the two dimensional Navier- 
Stokes equations can be modified by writing the continuity 
equation (3.18a) as: 

9u 


1 ^ 
r 9r 


(jO r- (rv) + 0) -r— = 0 

d Z 


and adding it to the left hand side of eq. (3.20). Then the 
so-called "conservative form" of the vorticity transport equation 
is obtained: 


9w , 9 
9t 9r 


(ojv) 
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(3.21) 


Introducing 


1 ^ 
r 9r 


the incompressible stream 
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i ii = _ 

r 9z 


V 


function i|j defined by 

(3.22) 


into the continuity equation (3.18a) gives the following stream 
function-vorticity relation 
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(3.23) 


Because the magnitudes of the stresses are often very important 
in the fluid field (they are, for example, responsible for 


17 


OF FOOH 

generating the turbulence), the components of the stress tensor 
for the incompressible axisymmetric flows studied here are of 
interest. These stresses in dimensionless form and in terms of 
velocity gradients and the laminar Reynolds number are: 


rr 


2 _ ^ 
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2 3u 
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(3.24) 


Recent results for the two-dimensional driven-cavity problem 
indicate that convergence was more rapid v;ith the formulation 

. . ri2i 

than with the primitive formulation that includes Uie pressure. 

This study found that the accuracy of the primitive variable 
solution was very sensitive to the convergence tolerance used 
in solving for the pressure. Therefore the t|;-w equations are 
used in this study to predict the flow field variables. 

By. substituting the vorticity given by eq. (3.23) into 
eq. (3.20), a single equation in terms of the stream function 
is obtained. This equation, called the "biharmonic equation", 


is : 
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3t 
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2 4 

where the operators .D and D are defined to be 


2 2 
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(3.26a) 


and 


4 2 2 

D^!p = D^(D^tjO 


(3.26b) 
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Although the biharmonic equation contains only one 
unknown ip , it is a fourth order nonlinear partial differential 
equation and is usually more difficult to solve. In [11] 
an efficient method of solving the steady planar biharmonic 
equation has been suggested. However, this solution method has 
poor convergence characteristics for a turbulent field, 
the two equations for the ip-w relations are considered 
hereafter. 

Another minor variation of eq. (3.21) has also been used 
together with eq. (3.23) . These forms differ only by the 
manner in which the diffusion term is v/ritten. For e.xample, 
if in eq. (3.21) the 

i i_ rr 

r 3r ^ 3r^ 


term is used instead of 


9 , 1 9u) 


then, the following equivalent form is obtained: 
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Really, there does not appear to be any significant advantage 
in using one of these equations rather than the other, since 
neither is truly in conservative form. The physical, inter- 
pretation of the conservative form of the equations for fluid 
dynamics has been discussed in the past.'' ' If, for example, 
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the total flux of vorticity is conserved for a finite volume in 

the region of interest, then the vorticity transport equation 

is said to be in conservative form. Equations (3.21) and (3.23) 

have previously been referred to, by inspection, as conservative 

forms by several investigators. This misconception results from 

the natural transfer of vast experience with the planar equations 

to the much less frequently used axisymmetric equations. In 

other words, the identical procedure of adding a modified version 

of the continuity equation to the convection terms of the 

vorticity dynamic equation is used in both planar and axisymmetric 

cases. In both cases, this puts the convective terms into the 

conservative form without affecting the diffusion terms. In the 

planar case, this simple manipulation produces the conservative 

form desired. In the axisymmetric case, however, this is not 

^ 9 (x) 

true since the term — tt— in the diffusion oart is not in the 

r dr 

[14 ] 

proper form. This point has only recently been discussed 
for curvilinear coordinate systems. 

3 . 6 The Mean Field Equations to be Solved 

In the light of this discussion tlie following form of the 
nondimensional system of equations will be considered in this 
study because they have been found to have the best conservation 
properties: 


|- (- 1^) + 1^) + rQ - 0 

8z r dz 3r r 3r 
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where the new variable Q is defined to be 


fi = 
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and Sjj is the source term of J2 defined to be created only from 
the turbulent variables: 
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The equations for the turbulence model’s quantities, k-e, are 
given by eqs. (3.5) -(3.8) . The conservative form of these 
equations is: 
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where G is the generation of turbulence given in (1.7) a? follows: 

G = (— + + 4[ (lX) + (Z) + Z ilZ] 

^ ^3r 3z^ ^^4r^ 'r' r 3r^ 


(2 If- Hi 


(3.33) 


. (21 

Suggested values for the various constants are given in 
eq. (3.9) . 

The final forms of the temperature and the species transport 
ecuations are: 


|- (ruT) + |- (rvT) 
oz or 


i_ (r\ 11) + i- (rx 11) 
3z 3z^ ^ 3r ^ ^ 3r^ 


G + 

ReE T 


(3.34) 


H + h = li lf> 


(3.35) 


where 


^ " RePr ■*■ Pr, 


(3.36) 


D = 


1 t 


ReSc Sc, 


(3.37) 


= -A [1 - I (2z - 1)^] (1-r) 


(3.38) 
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ncmeri::al prediction of tme flow field 


•i . 1 FinifcG Diffaraneo Approximations for the Derivatives 

Suppose that the finite two dimensional or axisymmetric 

region is subdivided by two families of lines parallel to the 

two coordinates z,r. For any point P » (s ,r ) / eight 

P P 

neighboring nodes are considered; they are located on the hori- 
zontal and vertical lines around P, not necessarily equally 
spaced. The eight neighboring nodes are given compass 
abbreviations, as seen in Fig. *1.1. Thus, SE stands for the 
point (Zp + hg, r^-hg), etc. Let 4>(NE), MP/f etc. denote the 
values of ''>(z,r) at the points NE, P, etc. 



♦ I 


E\ 


Hs 


. / 

>.S6 


Fig. 4.1. Mode Abbreviations for a Non-Constant 
Spaced Norroal Grid Mesh. 


An approximation to the first derivative at the point P can 
be obtained by considering the. Taylor’s series with a remainder 
in the two points E and W, 


'E 


V 




1 

6 


A 

inp rr 
^ .O 




Z < 

P 


< 


z 


E 


(4.1) 
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(4.2) 




'w “ “ ^Zp Nir ’*’ 2 '^zz N‘/ ~ 6 ’^zzz^*’w^ ^w 

P p 


z < ;; < z„ 

w w P 


Solving (4.1) and (4.2) for •{)„ eliminating (|),^ assuming that 

Zp ZZp 

is continuous, one can get: 

zz z 


w 


h„-h 

^ .{)„ - 


h. 


h„(hp+h ) "E ^ hph^ 
E E W E W 


*■ ‘'„(hg+h„) *„+0<h^) 


(4.3) 


where h = max(h., h„) . Equation (4.3) has a truncation error OF 

w b 

( z) 2 

the first z derivative R of the order h with the form of? 


(z)_ ^e\ . 

^zzz^'’^ 


"w 1 ^ ^ "e 


(4.4) 


The analogous formula for <{) is of similar form and obvious . 

P 

This formulation is from the centered differences for ct> , which 

z 

is far better than that from the forward differences which can 
be obtained from eq. (4.1) 



^E~*P 

^E 


(4.5) 


or from the backward differences which can be obtained from 
eq. (4.2) 


^w 

h 

w 


(4.6) 


These are only first order accurate. 

The analogous second derivative may be obtained by similar 
Taylor expansions, like eqs . (4.1), (4.2) carrying one term further. 
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4 . 2 Finite Difference Approximations to the Governing Equations 


A quantitative description of the flow field under consi- 
deration is obtained by replacing the various derivatives in the 
governing equations with finite difference approximations and 
solving the algebraic system of equations. Before doing this 
the method of solving this system should be determined. It does 
not mean that a choice between a direct solution method or an 
iteration method has to be made first, but rather that the way 
of coupling the different variables of the field should first 
be considered since it will affect the form of the finite 
difference approximation. The finite difference models of the 
advection in the governing equations will be discussed first 
followed by the method of solution. 

4.2.1 The Convective Terms 

It is known that central difference approximations of the 

[35] 

convective terms may give rise to instabilities. These can be 
eliminated by employing one-sided difference schemes which 
assure that the flow numerical information is consistent with 
the physical flow. This improvement was introduced in 1953 and 
is discussed by some authors.^ ' ' ^ This idea was incor- 
porated in the definition of the " transportive property" 

which states, "A finite difference formulation of a flow equation 
possesses the transportive property if the effect of a pertur- 
bation in a transport property is advected only in the direction 
of the velocity." The finite-difference expressions of the 
convective terms which satisfy this property are called upwind 
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or upstream differences and depend upon the direction of the 
velocity component at each point. Two methods for carrying out 
this idea were formulated. In the first method, backvrard 
differencing for the advection terms is used if all the 
velocities in the neighborhood of a typical point P are positive, 
and forward differencing if all the velocities at this point 
are negative, that is; 



r ^ 

f ‘“p 

- 

^P 

n 

u 

w 


if 

Up > 0 


1 

'• E 

■^E " 

n 

4>pl/hE 

if 

Up < 0 



n 


’ if 

if 


Vp i 0 

Vp < 0 


(4.12) 


A similar formulation may be obtained by using the control volume 
approach which provides a conservative and transportive differenc- 
ing method to handle these equations. Unfortunately, the increase 
in stability gamed by upstream differencing is at the cost of 
accuracy. The truncation error of the convective terms is 
increased^ to the first order, reducing the formal accuracy of 


the method to first order. Upstream differencing has gained 

popularity in spite of this critici.sm, and it has been found, 

from a practical point of view, to recover almost the same order 

of accuracy as central differencing. The upstream differencing 
[171 

was found (by comparing about six different explicit methods) 
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to be the best compromise between accuracy and computinq time. 

It was also foimd^^^^ that the upstream dil ferencing of con- 
servative forms yielded results which were as accurate as those 
from the second order difference equations. 

Another upwind differencing method is known as the "donar 

f 19 1 

cell” method. This technique is similar to the above upwind 

technique, with the added advantage of retaining the property 

of second-order-like accuracy of centered space derivatives. 

However, it can be shown that in axisymmetric form the accuracy is 

still reduced to the first order. The K-R technique for differonc- 

# 

ing the convective terms is known to be second-order accurate 


in the converging state, and unconditionally stable during the 
course of iterations. This is because it is an upwind 
differencing technique with a correction to recover the second 
order accuracy which is treated explicitly from the previous 
iteration (or time step) . For example, the z convection of 
the quantity ijs will be differenced as; 


(u^^) ^ = Y, 


3z 


n ,n 


h 


E 


n . n 
“p '"p 


+ a->,) 


n 

V 


u^ i{' 


n . n 
u 

w w 


h. 


w 


h 


Wv r.n-1 


" - If. h^> “ 


(4.13) 


where 


h^+h 


iT 

w w 


“w*w - “pOpl 

w 


(4.14) 


and 
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= 0 if Up ^ 0 

= 1 if Up < 0 

g 

similar expressions are used for • 5 — (vcfi) , with a parameter y 

0 XT IT 

instead of y in eqs. ( 4 . 13) - ( 4 . 15) . The upper indices n and n-1 
z 

refer to the appropriate time levels. In the steady state 
(convergence state) and eq. (4.13) represents a 

g 

finite difference for the first derivative - 5 — (uct)) of second 

0 Z 

order accuracy. 

4.2.2 Modeling of Non-Linearities and Coupling 

One of the main differences between the approach using the 
primitive variables and the approach using the vorticity-stream 
function variables to obtain numerical solutions of the two 
dimensional flow equations is in the treatment of the advection. 

In the primitive variables approach the convection is non-linear 
since products of the variables and their derivatives 
appear. In the regular (or ip-fi) approach the velocities appear 

only as coefficients of the vorticity derivatives and the ip and w 
variables appear in a sort of quasi-linear form. 

As has been stated in many works the 'Jj-co solution method 

has to be under-relaxed in order for it to converge (also with 
upwind differencing) and this is due to the instabilities that 
the boundary condition for the vorticity introduce into the field. 
Since the boundary conditions for to consist of linear relations 
between ({; and to [see Section 4.3.2] it is reasonable to assume 
that the coupled solution of the two equations for ip and to will 
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result in suppression and possibly elimination of the instabilities 
A simple way of coupling the equations is to leave their 
finite difference forms as in the separate case, where the 
coupling between tp and w is done only on the boundaries . A more 
intensive way of doing it is by coupling the and w variables at 
every point of the field. '■ In this case the following second 

order approximation for the nonlinearities of the convection in 
the non-conservation form is considered: 


r 9 / \ 1 n 

lu ^ (o:)] 




(4.16) 

where n is again the time (or iteration) level. Applying the 
continuity equation (3.1) and the stream function definition 
(3.22), the following equation can be obtained for the planar 


vorticity convection. 


[9] 


r„ . 9oJ, 


n 


I- (u^-^ to’^) + I- 


II 


n-1 


] - 




III 




IV 


V 


VI 


(4.17) 


In the above equation, terms III and VI are source terms known 
explicitly from the previous time step. Terms I and IV represent 
the convection of the vorticity and may be treated by the K-R 
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method presented in eq. (4.13). Terms II and IV are nev? 
terms which represent the convection of the stream function in 
the vorticity equation. It can be easily shown that the 
diagonal dominance of the coupled system as well as the stability 
conditions will not be affected if the terms II and V will be 
central differenced or upwind differenced. 


4.2.3 The Governing Equations 

As was mentioned in section 3.6, the quantity fj defined by 



has "better” conservation properties than the vorticity oj 
itself. The relation is given in eq. (3.28) and the dynamic 

equation for is given by eq. (3.2a) , or by the following 
equation in the non-conservative form 


+ ''If = - 


3z' 


i]xa) + ^ im + 1 (ys^) + 




(4.18) 


8r 


where y = + y^ is the total viscosity coefficient. • Let 

us define the following grid parameters at the point P: 




"r = 2 ' 


a 

z 


hj,/h 


/ 


a 

r 


h^/h 


S 


(4.19) 


The finite difference form of the stream-function vorticity 
relation, eq. (3.28), will be 
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(4.20) 


This equation is taken alv/ays to be at the new time level, 
say, n. 

The finite difference form of the vorticity dynamic equation 
(3,29) will be: 

* , 1 
*Nt2h^ 
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This equation is solved at the n time level where all the coeffi- 
cients u, y and the source terms , Dj- and R are taken in the 


n-1 time level. 
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4.2.4 Tne k-e Sanations 

The solution of the k-e equations, (3.3i~3.32), is obtained 
in this study in a coupled manner, based on a known flow field 
mean velocity distribution. The nonlinear terms are quasi- 
linearized with respect to the time level (or iteration index) n: 
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(4.25a) 


(4.25b) 
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(4.25c) 


The finite difference approximation of the k equation (3.31) 
is : 
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(continued) 
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and the finite difference approximation of the e equation (3.32) 
is ; 
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where : 
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= (U,H. jij/p 


V p2 (n-l) 
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(4.27) 

(4.28a) 

(4.28b) 

(4.29a) 

(4,29b) 


The operators D and D are defined in eqs. (4.23) , and the 

Z 2T 

generation term, G, is defined in eq. (3.33), taking the velocities 
from the last time level,- n-l. 
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The finite difference approxiirations for the temperature 
and species equations (3.35~3.36) may be obtained in a manner 
similar to that used for the above equations. The equations are 
linear and may be solved in an uncoupled manner. 

4 . 3 The Boundary Conditions 

The system of equations governing the flow field, {3.28-- 
3.38), are of an elliptic nature and therefore boundary 
conditions for all the variables should be specified on all the 
boundaries. Four boundaries confine the present field: the 

outer cylinder wall, the axis of symmetry and the inlet and 
outlet cross-sectional planes. Although the main purpose of 
this study is to solve the turbulent flow, the program was 
checked with some laminar test cases and therefore the laminar 
boundary conditions will be also mentioned. Let us begin with 
the symmetry axis conditions. 

4.3.1 The Axis of Symmetry 

Along the symmetry axis the stream function has a constant 
value and it xfi convenient to choose it as zero: 

•Mr = 0) = 0 (4.30) 

so that ijj simulates the mass flux contained between r and 0. 

The values of 9, can be obtained from its dynamic equation (3.29) 
assuming symmetry, i.e.: 

V = 0 (4.31) 

on this line. Therefore in the limit as r 0 that equation 
will assume the following form: 
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(4.32) 


(4.33) 


It ifi important to note that other formulae^ ^ for 
can be obtained bv' annuininq abov.? eq. (4.3)) thot t’lo g darivativen 
in ocr. (4.32) are no<i'' idiblo, and tlu5re.l.or(' with = 0 
wo can qot that 


s)(r 0) = const - C 


(4.34) 


whore the value of C can bo evaluated from the il)-d equation 
(3.28) with 


(r " 0) = 0 


(4.35) 


since the limit for the centerline u velocity definition 


, . 1 3 il' ^ 

u = lim - 

r-^-O ^ ^ 3r r=0 


(4.36) 


exists. iM the pia^yont approach, both equations for 

’!» and d are solved exactly on the axis of symmeti-y. Since 


„ 3C .. ;^k _ :)i 
i)r ” Or 3r 3r 


(4.37) 
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on the axis of symmetry the dynamic equations for k, e, T and C 
on this line may be derived from equations (3.31), (3.32), 

(3. 34), (3.35) giving the foliow..ng set: 
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(4.38c) 


h = h <° H> + ° 7^ 

d Z 


(4. 38d) 




(4. 38e) 


4.3.2 Solid ?7all - Laminar Case 

The stream fimction is constant on tl\e wall since the mass 
flux of the flow does not change during its passage through the 
system. Concerning the laminar vorticity on the walls, with the 
classic approach the dynamic equation for n is not solved on the 
rigid wall surface, but rather the values of are taken from the tji-fi 
relation (3.28). By using the no-slip conditions, (|^) • = (~^) =0, 

rr Y, 

w w 

the following second order accurate values for may be obtained 

. (221 
by expanding i|i and fl xn a Taylor series about the wall point w: ^ 


(4.39) 
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Fig. 4.2: Grid Near a Cylindrical V’all 


and, for z = const .is shown in Fin. 4.3 
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Fig. 4.3: Grid Hoar a Vortical Wall 
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W 1 

h r 
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(4.40) 


where 1 is one point away from the wall point, w, into the field and 
h is the distance between the points 1 and w. 
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A different approach to get Q on a r = const. = r wall is to 

w w 

impose the no slip condition (4.39) in the relations to get 


n = . 1 (i!t, 

r^ 3r w 


with 


(— ) 


= 0 


w 


v/hich gives; 


“w ' -2^ <*1- V 

W 


(4.41) 


For a z = const, wall one can use the no slip condition to get 

(H) =0 

w 

And since ip is constant along this wall an equation the same 
as (4.41) holds true. 

4.3.3 Solid Wall - Turbulent Case 

The main difficulty in imposing the turbulent boundary condi- 
tions on an impermeable wall is the matching between the 
turbulence model and the boundary condition formulation. Let 
us first discuss generally the nature of turbulence near a wall. 

The above recommended k-e model was developed based on certain 
assumptions concerning isotropy of e and the relations between the 
Reynolds stress and the turbulent kinetic energy k, among 
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other assumptions. Those assumptions are not valid near a v/all 
which has the following contradictory features: 


(i) sharp flow properties variation. 

(ii) the turbulent viscosity begins to affect the various 
turbulent processes due to the low turbulent energy 
level (or low turbulent Reynolds numbers) . 

(iii) the presence of the wall destroys any isotropy 

feature that is carried to it from the main stream. 
There are two commonly employed ways of accounting for the 
turbulent phenomena in the wall region: (i) the v/all function 
method and (ii) the method of modelling low turbulent Reynolds 
number phenomena. It is very common to divide the turbulent 
boundary layer with thickness 5, into two main parts: the 

V V 

inner layer, 0 £ ^ £ 0*2, and the outer layer, 0.2 £ j £ 1.0. 

The inner layer is dominated by the wall effects and the outer 
layer by the flov/ effects. Since the k-s model was designed 
to account for the outer layer properties, the inner layer 
determines the turbulent boundary conditions. The influence 
of the inner layer depends on the relative size of the terms on 
the right hand side of the wall shear stress 

where u‘ , v' are the turbulent fluctuations in the longitudinal 
and lateral directions, and - p u'v' is the non-negligible 
Reynolds stress in the boundary layer. It is very convenient to 
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V t 


yF poc>R u / 


4* •{• 

use the U and y formulation in the inner layer as follows: 


V* 


y V 

V 


(4.43a,b) 


where v is the shear velocity defined as 


V = /t“/p 


(4.43c) 


Three regions may be distinguished in the inner layer: 

(i) The laminar sublayer: 0 < y"^ < 5. This is the 

' innermost part where 

3U 


- u ' v' << y 


9y ' 


and 


(4.44a) 


+ + 

U = y ; 


(4.44b) 


(ii) The logarithmic law region: 40 < y”*" < 0.26'*'. 

This is the outermost part where the only contribution 


to the stress is -p u'v', and 


4" 1 4" 

ir = i £n y^ + C , 


< = 0.41 and C ~ 5.1 


= ^ Jin (Ey'*') , E :: 9.0 25 


(4.45) 


(iii) The buffer region: 5 < y <40. This is the inter- 

mediate part between the sublayer region and the 
logarithmic region. 

In the "wall function" approach, the first grid point is 

located away from the wall, and the goal is to locate it in the 

. . + 

Ipgarithmic region. In this region and k are constants, and 


k i>v jui |f3^ 

OF Pi>0/C 0»?/iXE7Y 

the generation balances the dissipation in the k equation (3.31): 


V 


_ _l/4 .1/2 


c3/4 j^3/2 
- ^ — 


(4.46a) 

(4.46b) 


U. = /k K Y 

where y is measured from the wall, and k is the Von-Karman coeffi- 
cient ( ~ 0 . 41) . 

From the balance between the e generation and dissipation in 
eq. (3.32) one can get: 

.2 


^2 - ^1 = 


K 


r cl/2 

e y 


(for y >> 1) 


(4.47) 


This equation is the only condition which assures that the k-e 

model also holds in the logarithmic zone. Generally and C 2 

as well as are constants found from the experimental data of 

F21 

one-dimensional turbulent grid flows. It can be seen that the 

constants appearing in eq. (3.9) do not fulfill eq. (4.47). 

In [7] another set is recommended: 




= 1.47 ; C 2 = 0.18 ; = 1 ; 


Cg. = 1.3 (4.48) 


This set of constants is not in agreement with the one that was 
suggested for obtaining a good approximation to e in the main stream 
field. Analyzing the turbulent quantities in the sublayer region, 
one might conclude that 

k**" ~ (y^) ^ = constant . (4.49) 

Eqs. (4.49) leads to the following relationship between the model 
coefficients 


^ (2C3_/C C 2 ) = y 


(4.50) 
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suggiist a dependence of the model coef- 


Therefore some papers 


ficients on the local turbulent Reynolds number defined by 


^t = 


C k' 


(4.51) 


or on the y normal distance. Such a variation can be obtained by 


adjusting the constants using some experimental data: 


. [ 1 ] 




Cj = [1 “ 0.3 exp(- 0.125 • y )] 


125 


^]i 50 + rJ 

t: 


(4.52a) 

(4.52b) 


It can be seen that these variations doe not resolve the above 
conflict since in the logarithmic zone the Reynolds number is 
highly turbulent. The above functions do allow us to locate the 
first grid point in the sublayer region, however. 

The following set of constants was adopted in the present 
study to get good agreement with the axisymmetric experimental 
results : 

= 1.5 ; C2„ = 2.0 ; = 0.09 ; = 1.12 (4,53) 

With the "wall function" approach the first grid point is 
not located on the wall but at a distance y away from the wall, 
as is illustrated in Fig. 4.4. 



OP 


V-' f--l 

vii/ji.jTV' 





Fig. 4.4: Grid Near the Wall for a Turbulent Flow. 
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For this point p, located in the logarithmic region v;here 
is the sublayer edge, the following relations exist: 


+ V 

y = y - 


u+ = , U-" = i Jn(E y+) 


(4.54a,b,c) 


= I ^o " I ^o - 1] [i^n(E y-^) - 11 


~ y'’' (u"*" - 2.5) - 38.5 


U.+ = „ _1_ = . i.1 

tr + + 

Ky y 


(4.54d) 

(4.54e) 


K = 


e = 


k’*' = 3.15 
o 

3/4 +3/2 

^la ^ 1 

K H 


2.24 


(4.54 f) 
(4.54g) 


E = 9.0 


(4.54h) 


If the first grid point happens to be in the viscous sublayer, 
then: 


+ 

u = 





(4. 55a,b,c) 


K 


0 


1 


+2 

y 


= 0.2 


(4.55d,e) 


If the first grid point cannot be kept in the sublayer or 
logaritlimic regions during the iterations, the approach 
of continuous v/all functions can be adopted. In this case 
the first grid point may be located on the wall and the following 
boundary conditions have to be considered for the first grid 
point away from the wall: 
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' 

OJ- ip '.£',Y 


- xR 


A y'*' u"^ + (1 - 0 ik y’^(u'*'~ 5) - 3n.v.] 


= -tl + (1- e - 1) ] 

y 

4. 4-^ 4.2 

k = 0.1 i' + (1 - G (3.15 - 0.1 i' ) 


C = 0.2 + (1 - 


-uR 
G ( 


+• T IX 

X. ^ ^ 4U • 4 ^ U 

+ 


» 0 . 2 ) 


(•1.56a) 

(4.56b) 

(4.56c) 

(4.56a) 


wliere 


a “ 2.5 


(4.56g) 


x^ras chosen to match the experimental u orofiles 


r ^ 4 1 

In recent years some carefully measurements of scalar 
quantities in the wall region have been carried out. The profiles 
of s'** in the logarithmic x*egion was found to satisfy the 
correlation 

= ~~ ?n y*^ + F (4.57a) 

with 

Kg =0.46 

where s’** is the dimensionless scalar (temperature or species) 
defined to be 


+ 

S 


S- S 




w 


(4.57b) 
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For the case of the species S = C/ the boundary condition is 

(|f) = 0 (4 

^ w 

and the Schmidt number is taken to be 


Sc^ = i- = 0,885 (4 

^ ^c 

The turbulent Prandtl number was found to vary through the 

[25] 

boundary layer region. The first suggestion which was 

f 2 G 2 *7 1 

based on some measurements ^ ^ was : 

« 

2 

Pr. = 0.9 - 0.4 (^) (4. 

t 0 

r 2 8 1 

and was updated later ^ ^ to 

2 

Pr^ = 0.95 - 0.45 (p (4 


Since the adiabatic wall is not of interest in this study, the 
following temperature continuous wall function will be adopted 


T- T 
w 


+ 


-aR 


/X (^) 


= y + (1 - e 


[2.15 £n(12 y'*') - y"^] 


(4. 


w 


Pr 


t 


0.95 - 0.4 (p 


(4. 


and for the species equations (4.58a,b). 




58a) 


8b) 


59a) 


59b) 


60a) 


60b) 
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4.3,4 Inlet Conditions 

In most problems the flow velocity temperature and the 
various concentrations are known in the inlet section i. Usually 
a fully developed velocity profile is given from which and Q 
may be deduced: 
r 

’J' = / y U. (y) dy (4.61a) 

0 ^ 

(4.61b) 


For the laminar case u. is taken as 

X 


u. (r) = (1 

X o 




(4.62a) 


and for the turbulent case the velocity defect law is adopted: 
for highly turbulent Reynolds numbers: 

.+ 


+ + 1 r' "*^^t + ""^^t 

- u+ = (| An 5_) ( 1 - e ^) + y^ e 

y 


(4.52b) 


for low turbulent Reynolds numbers: 
“o 


(4.62c) 


For the turbulent energy and the turbulent dissipation, we 
assiome that v/e have the profiles that are described in equations 
(4 .56c,d) . 


OF pio^ 


In order to have a close form of the system at the inlet, 

we have to specify the ratio 

* 


6 


U 


(4.63) 


For laminar flow it can be shown that 


3 = / 


Re 


(4.64) 


and we shall assume similar 0 for the turbulent flow. In any 
event the inlet conditions have very little effect on the inner 
field solution, especially if the inlet section is taken to be 
far from the regions of high gradients. 


4.3.5 Exit Conditions 

At the outlet boundary we assume that the longitudinal 
diffusion is negligible, that is 


2 

8 r all the y 

- 2 ^variables’' 
9z 


(4.65) 


and the flov^ is basically parabolic near the exit section. 


4.3.6 Inner Tube Trailing Edge Conditions 

The inner cylinder is indxcated by the point p in Fig, 4.5. 

i 9''-r--p H- 



Fig. 4.5: End Wall Computational Region. 
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One approach ' that is easy to use and most coirunon to consider 
that the dividing stream line, has no curvature between 
point p and point E, satisfying, at the end wall. 



The approach that will be used in this study is based on the 

physical fact that has two different values at the point p, 

depending on whether the vr derivative is calculated from values 

above the wall (to give or from points below the wall (to 

give . Those two values are opposite in sign. Solving the 

n equation at the point N, v/ill be used. Solving at the 

P 

point S, the v;ill be used. When solving the equation at 
P 

the point E, n = o will be used since ^ is a separation line, 

P P 
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The solution to the set of six finite difference equations 
is obtained in this study in an iterative manner since it is 
impossible to invert the entire system directly. We use the 
term, one global iteration of the flow field, here to mean that 
the following tv/o piece chain is completed: 

(i) For a given k-e field the mean field quantities, 
i.e. ip~Q, are obtained to a certain convergence 
tolerance, including the velocity fields. 

(ii) For these velocity fields the k-e -turbulent variables 
are obtained iteratively to a certain tolerance. 

After the global iteration set converges, it is possible to 
solve the concentration field, C, and after convergence of 
this field it is possible to solve the temperature field. This 
iteration process may be sketched as follows: 
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5. 1 Solution Technique 

The solutions of the scalar fields C and T are obtained by 
the successive line relaxation (SLR) technique. With this 
technique the scalar variable is solved implicitly for all the 

rows and then for all the columns in the field, replacing 
successively the old values by the current nev/ values. The 
solution of the coupled systems, like il/-n and k-£, are obtained 
by the successive block line relaxation (SSLR) , which is the same 
as the SLR technique but expressed in terms of blocks instead 
of scalars. 


5 . 2 Stability Analysis 

The basic blocked system that is solved here may be described 
by the following equation 

[E] + [w] + [N] + [S] + IP] + R = 0. (5.1) 

where [E] , [w] , [N] , [S] and [P] are the coefficient matrices 

(in the present study 2x2), R is the source term vector and 
^ is the variable vector. The SBLR along the field columns may 
be written as 

E^"^ + w”^"^ (|)JJ + <j)g + Pj"^ 

+ (R^"^ + p2"^ = 0 (5.2) 

where it is understood that E is a matrix and P is a vector. 
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The coefficient matrix [P] of the middle point is split 
into two matrices 

[P] = [Pj] + [P 2 ] 

where [P^^] contains the total diffusion effect and the convection 
in the implicit direction [which is r in aq. (5.2)] and [P 2 ] 
contains the effect of the convection in the other direction. 
Treating the [P^^] contribution implicitly and the [P 2 ] contribution 
explicitly has been found to be the best way of splitting of [P] 
for turbulent flow. For laminar flow, it is found to be best to 
set [P 2 ] = 0. This is due to the special boundary condition near 
the wall used with turbulent flows. The stability requirement 
is expressed in the diagonal dom^inance inequality: 


, |w| + 1 e| 

I Pj^ I and (5.3) 

^ |S| + |N| 


where strict inequality has to be maintained at least in one 
point of the field. 


5.2.1 Stability of the System 

Let us check first the stability in the z-direction of the 
system, eq. (4.20-4.21): 


E + w = 


r ( 1+0 ) 
= JE- z 


4 h^o 
z z 


v„ — + V + u h 
E 0 w w z 

z 

^E I^e'^z 

+ V + = 

z z 


if u > 0 


if u < 0 j 
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+ 1 + ^ 


z z 




"'p 1 , '^■^°z>^ 1 3 

2 a, ^2-2h^rp 'r 


(o.- 


l+a 


u 


ji r 

2K p 


p 


if u > 0 
if u < 0 


r A 
* ^ 

r 


B 


— ) + ^ (a 
a h z 

r z 


-) ] 


so that stability requires: 


if u > 0 


+ u 

w z 


if u < 0 


(y^ z 


-±1 + V 
a w 

z 


1+0 h 2 

" <ir> ' 


1+a. 


- + 


2 / T ^ Nv 1 , S 

a (1 + —) 1 + — 

r rp Tp 


“p \ 


if u > 0 


ie u < 0 


h 2 1+a 

^ 4(l^a^) (j^) 4- - 

IT 17 ^ ^ 


a (1 + ;r 

r r, 

^r 1 

-^] Vp + (- + 1) Vp 

o Z 


1 + 


QIH 





It is reasonable to assume that o and a are analytic 

r z 

1 2 

f\inctions of r and z such that a - — = o(h ) and therefore may 
be neglected in the stability analysis. It can be seen from 
eq. (5.4) that this scheme is unconditionally stable and it is 

[91 

more stable than in the case where '\i and R are solved separately. 

It can be shown, also, by checking the spectral radius of the 

system, that the rate of convergence of the coupled system 

* 

R_ is greater than that of the separate system R_ : 


R. 


= 1 + 


c {l+a ) , 

■ rcuij (Re^ + D) 

i Z (J ^ 


- + 


(h^) 


(5.5) 


where Re^ is the averaged cell Reynolds number 


Re^ = (Up h^ + Vp h^) Re (5.6) 

3 Q 3 n 

D is a Dositive function of the 9, derivatives -:r— and and 

*■ 3 z 3r 

2 

is 0(h ). Therefore, it may be concluded that the main 
improvement in the rate of convergence (and also in 
stability) of the coupled system results from treating the 9 
in the ip-Q relations as an unknown rather than as a source term. 


* Rg was taken to be the rate of convergence in the case that 
for every 'p solution, an fi solution is obtained. In that 
case it may be proved that r| ~ R^j^ • R^ where R,j^ and 

are the rate of convergence of the equation and of the 
equation respectively. 
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wva _ j. 1 ^. 

Ti'CvUtinvj the convection l?y aplittinc it between •;• and d ai\vi in 
eq. C4.17) has second order (h“) effects on the rate of convercjence . 
But since it is hard to evaluate the coefficient of this 
effect, this split is used to obtain the results of this effect. 

No extensive study of this effect has been done but by using 
this technique a higher rate of convercrence can be achieved. 
Generally, computer results have shown that the actual rate of 
convergence is higher th.an is predicted by eq. (5.5), since the 
boundary conditions are based only on the c-d relations. 

5.1.1 Stability of the k-t System 

The stv^bility condition for the k-t syste.m can be derived 
in a manner similar to that of the 'v'-d system. It turns out 
that the upstream influence ma.kes the system more stable. Vie 
consider first the 'aero velocity case. In this case the 
matrices E and W are po.sitive definate, and therefore if 
is also positive definate we liav’e: 

4 h*l - w w p p Iw^^ 


s 


'"p i f "p) '=3 ^ =1 = 


h - 


e h‘ 


- C,j(4C2-Gi) G (5.7) 

It may result that the oossible critical points in the field 
for which this conditions should be chocked, are the points in the 
logarithmic 'oone. Substituting the locarithmic relations, taking 
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into account only the leading terms (of y powers) , the following 
equations can be derived in terms of the "+" system, with 


^3/4 j^+1/2 

^ > 4 = 2.85 


(5.8) 


where R is the outer cylinder radius and h^ is the distance of 
the second point away from the wall. This condition always holds, 
since we usually have more than four points in the r direction. 

If the convection is not negligible, then the condition (5.8) 
would be less severe. Taking into account the convection in 
the logarithmic region it can be shown that 


2.85 + £n(Eh'*') 

1 + S,n(E h"^) 
z 


(5.9) 


The ratio of the rate of convergence of the coupled to the non- 
coupled k-e system can be shovm to be approximately, for 


2.85 + ?.n(E h'*’) 

~ 1 + 5 

^ + 21.8 £n(E hp + 9 

which is about 1.15 r 1.18 for practical solutions. 


(5.10) 


5 . 3 Application of the Turbulent Wall Boundary Conditions 
The wall boundary conditions for the turbulent field 
solutions are not applied implicitly in the computational 
procedure. After every iteration the near-wall values are 
updated as follows: 


59 


OF POOR QUALtTY 

At point 1, which is a distance h away from the wall 
(see Fig. 4.2 or 4.4), the new values of all the field quantities 
are known after the nth iteration. Then the following steps 
are taken: 

(i) solution of the non-linear equation for v*/Re from 
eq . (4.53a), 

(ii) obtaining oj'^(y=0) from (4.53b) and replacing the 

old n on the wall by the new w/r , 

w 

“f" *4" • 

(iii) solving eq. (4.53c,d) for k and e and replacing 

e , e, and k, 
w 1 1 

Similar steps are taken for evaluating the temperature gradients 

on the wall for every iteration of the energy equation. The 

above explicit treatment of the turbulent boundary conditions 

may be shown to be stable by plotting on a curve of (-w ) vs. 

some of these explicit iterations as in Fig. 6.9. Since this 

curve exhibits convergence (the limit " exists) , and since 

+ 

the system of equations demands increasing values for 
increasing w values, this procedure is stable. 
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6. COMPUTATIONAL RESULTS AND COMPARISON WITH PREVIOUS WORK 


6 . 1 Grid Development 

A variable spacing grid was chosen to allow for more grid 
points in the regions of large gradients which are located 
as part of the solution process in order to obtain an accurate 
solution with a small total number of points. Since the large 
gradients of velocity, as well as those of turbulence quantities 
occur near the solid walls, more points will be located in the 
wall regions. An analytic distribution is preferred in order 
to minimize the effect of the (a - coefficient in the finite 
difference forms in Chapter 4 . Since we are mostly interested 
in turbulent flows, a geometric variation of the grid spacing 
might be a good choice since the mean flow velocities are almost 
linear on the geometric transformation coordinates. 

Figure 6.1 illustrates the grid in the r direction. The 
inner cylinder is located at j = N^^ where the total number of 
points in the r direction are N. The radius of the point Q 
has the value of ^ (l+r 2 ) , with the integer point j. It is 
assumed that 

h = a 0 < k, < N- ™ 1 for D < N^ 

r ^1 —1—1 -^—1 

k 

h^ = a q 2 ^ 0 £ k 2 £ N 2 - 1 for £ j £ N^ + N 2 - 1 

(see Fig. 6.1) . The points N^^^ + 1 and N^ - 1 both are the same 
distance from the inner cylinder wall. The geometrical coefficients 


61 




r 1 
• } 





and q 2 should fulfill the following relations 


N 

^1 " q^ ^'^l " 


a ^2 

= 2 5 ^ (g/-l) 


( 6 . 2 ) 


For a given "a", equations (6.2) will give q^ and q 2 * Usually 
"a" is determined beforehand, i.e., to give the number of 
points we desire to be located in the boundary layer region, 
or to guarantee that this point is in the logarithmic 
region. 

Although this non-constant spacing grid is an attempt to 
raise the efficiency of the calculation (by using a smaller 
number of grid points) it is still far from the optimal grid. 
This is because in the core region far downstream from the 
inner cylinder trailing edge, there are many more grid points 
than are necessary to trace the fairly smooth variation of the 
variables. This is however, possibly the best transformation 
that can be applied to this field without using much more 
complicated curvilinear coordinate transformations obtained from 
differential- mappings . A sample of the coordinate system 
around the trailing edge of the inner cylinder is given in 
Fig. 6.2. Here q^ = 1.38 and q 2 = 1.3 for r^^ = 0.362, and the 
geometric series factor in the z direction is 1.8. 



6 . 2 The Laminar Field 


rj 


A program for obtaining solutions for the mixing of a 

confined coaxial flow was developed, based on the numerical 

methods presented. The validity of the numerical procedure 

was checked by comparison with other results in the literature. 

It was found that the present procedure yields results for an 

entrance region flow with an initial plug profile entering a 

[29] 

pipe that varies by O.l^S from analytical results. ‘ 

« 

It is convenient to define the ratio of the mean velocity 
of the outer jet to the mean velocity of the inner jet as a 
parameter 



*2 ~ fl 

1 - r^ 


^1 


(6.3) 


as well as the mean Reynolds number 


Re„ = ({).)• Re (6.4) 

m i 

Another Check was made for confined jet mixing using the same 
field parameters as others; 

-^ = 1.17 ; Re^ = 496 ; r, = 0.5 

The maximum discrepancy from the analytical results was about 1% . 
Here a 41 x 31 grid was used with ^ ^2 ~ 

second order accuracy of the scheme was verified with 5,^ = 6 
chosen as the outlet section to impose the parabolic conditions. 


✓ 
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Predictions of the intensity of the reverse flow, 

rnxxi rn 

for rj^ * 0.5 and ~ described in Fig. 6.3 for 

fully developed inlet profiles. There is not much difference 

between the = 3 and the z^ = 5 results. The results v/ere 

not stable above Re„ = 3200. In the range of 3200 < Re„ < 4200 

m ^ m 

the error in the iteration procedure did not grow, but it did 

not decrease either. Above a mean Reynolds number of "SOOO 

the iteration procedure blows up. The zone of separation is 

a function of Re^ and U2/Uj^. In Fig. 6.4 the locations of the 

separation and the reattachment points as functions of Re^^ 

are described. These results give the impression that perhaps 

there is an asymptotic bubble length as Re^ tends to infinity. 

For r, =0.5 and Re = 2000, the laminar stream function is 
Ira 

described in Fig. 6.5 for two cases: ~ ^2^^^1 ~ 0*02. 

In Fig. 6.6, the temperature field is described. The 
temperature difference between the inner and outer jet was 
the reference temperature. The temperatures rise on the axis 
of symmetry is of the order of A/10 in eq. (3.13) . The 
maximum rise in the temperature is near the end of flov; field 
where the region of maxim’um heat generation is located. 

6 . 3 The Turbulent Field 

Calculations of the turbulent flow fields were carried out 

using experimental data for the inlet flow profiles. Unlike 
[31] 

other work the vorticity source term (3.30a) is not 

2 

neglected since for 0*1/ ^ ^ is important in the high 

shear domains, especially near the dividing streamline and near 
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the walls in the viscous sub- layer zone where 7^v x (n is 

the normal to the shear layer direction) . Similar to other 
works ' the present study was done with r^ =0.5, 

Re = 72,000, and a pipe length, 1,^ - 20. Predictions of the 
stream function and the turbulent kinetic energy for ~ 

are described in Fig. 6.7. The temperature field of the present 
problem is shown in Fig. 6.8 for the constant A = 30. The 
effect of turbulent diffusion can be seen by the over-heating 
near the axis of symmetry at z ~ 5. That is due to the relative 
high v^. Due to a lack of experimental results for this field, 
the only check of the calculated results that was made is the 
global conservation of the variables. Mass was conserved 
to a 0.1% tolerance, while momentum and mean energy was conserved 
to a tolerance of 0.5%. The fact that the turbulent boundary 
conditions are, essentially, treated explicitly, reduces the 
rate of convergence. In Fig. 6.9, some of these explicit 
iterations are shown for a wall point at z = 3. Figure 6. '10 
demonstrates the difference in the convergence rate between a 
field where the point ahead of the outer cylindrical wall is 
located in the viscous sublayer, and a field v/here this point 
is located in the logarithmic region. 
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21x81 GRID 



(T) point 2 in the 
sublayer region 

@ point 2 in the 
log region 


00 

ITERATIONS 


FIGURE 6.10-- VARIATION OF THE ERROR 
IN VS NUMBER OF ITERATDNS 



7. CONCLUSIONS AND SUMMARY 


This work is an attempt to solve numerically for the field 

of a two dimensional turbulent flow using a "two equation" type 

model for the turbulence. Since the flow field considered in 

this study is the incompressible confined turbulent mixing of 

two coaxial jets with internal heat generation, the mean stream 

function and mean vorticity were chosen as the flow variables. 

The turbulence effects v/ere represented with the k-e model. 

A conservative coupled-variable finite-difference scheme was 

employed. The finite-difference algebraic equations were 

solved iteratively by successive block line relaxation. 

Conditions for obtaining stable turbulent solutions were 

formulated. The double value of the vorticity and the poorly 

defined variation of the turbulent variables at a sharp corner 

were treated through the generation terms. 

Although the coupling between equations and k-e 

equations tend to increase the convergence of the solution, 

the explicit form of the boundary conditions tend to decrease 

the rate of convergence compared to that of the laminar field. 

Since the gradients of all the variables in the logarithmic 

region near a wall is much smaller than they are in the sublayer 

region, the rate of convergence there is much faster. The 

2 2 

downstream face boundary condition, 9 /Sz = 0, is the correct 
one to be applied for the so-called "fully-developed turbulent 
flow" condition since applying a 9/9 z = 0 condition contradicts 
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the boundary conditions on the wall. Applying such a boundary 
condition to a turbulent flow field can cause erroneously 
large turbulent diffusion coefficients in the field. There is 
no effect of the inlet conditions if they are applied far 
enough upstream from the trailing edge of the inner 
cylinder (~ 3 to 5 R) . 
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NOMENCLATURE 


D 

D 

G 

k 

L 

P 

Pr 

Q 

r 

R 

Re 


Re 

m 

R^,R^ 


S 

Sc 

S 


T 


T 


u 

V 

z 


concentration function 
k-e model coefficients 

partial differential operator defined in eq. (3.26) 

mass diffusion coefficient 

turbulent energy production 

turbulent energy 

typical axial length 

pressure 

temperature Prandtl ntimber 
mass flux 

radial direction and coordinate 
radius 

rate of convergence 

radiation number 

Reynolds mmber 

mean Reynolds number 

truncation error 

source term 

Schmidt number 

heat generation rate 

temperature 

axial velocity 

radial velocity 

axial direction and coordinate 
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Y 

e 

X 


u 

^eff 


V 

a 

o 

T 

' 1 ^ 


k 

e 


\jj 




the thermal radiative diffusivity 

integer to choose the right differencing direction 

turbulent energy dissipation 

thermal conductivity 

laminar viscosity coefficient 

'turbulent viscosity coefficient 

the effective viscosity coefficient 

kinematic viscosity coefficient 

turbulent energy 'Prandtl number 

dissipation Prandtl number 

stress tensor 

stream function 

maximum stream function 

vorticity 

w/r 


Indexes 



inner cylinder 
outer cylinder 
effective 
tensorial indices 
laminar 
radiation 


t 


turbulent 


